#include <cmath>
#include <limits>
#include <iostream>
//#include <iostream>
class Function
{
public:
    double operator()(double _x);
    double diff(double _x);//导数
};

class EquationSolver
{
public:
    virtual double solve() = 0;//虚类，不可以实体化
};

//double func(double _x);全局函数，需再定义头文件
class Bisection : public EquationSolver//继承EquationSolver，由于虚类，必须使用solve，否则错误//二分法
{
public:
    Bisection(double _a, double _b, double _delta, int _M, Function _B_f);
    double solve();
private:
    double a;
    double b;
    double delta;
    Function B_f;
    int M;
};

Bisection :: Bisection(double _a, double _b, double _delta, int _M, Function _B_f)
{
    a = _a; b = _b; delta = _delta; M = _M; B_f = _B_f;
}
double Bisection::solve()
{
    double h = b-a;
    double u = B_f(a);
    double c;
    double w;
    int k = 1;
    for(k = 1; k <= M ; k++)
    {
        h = h/2;
        c = a+h;
        w = B_f(c);
        if (fabs(h) < delta || fabs(w)<std::numeric_limits<double>::epsilon())
        {
            break;
        }
        else if(w*u>=0)
        {
            a = c;
        }
    }
    // std::cout << "h = " << h << std::endl;
    // std::cout << "k = " << k << std::endl;
    // std::cout << "c = " << c << std::endl; 
    return c;
};

class Newton : public EquationSolver//牛顿迭代法
{
private:
    double x0;
    int M;
    Function N_f;
public:
    Newton(double _x0, int _M, Function _N_f);
    double solve();
};

Newton :: Newton(double _x0, int _M, Function _N_f)
{
    x0 = _x0; M = _M; N_f = _N_f;
}

double Newton :: solve()
{
    double x = x0;
    double u;
    int k = 0;
    for( k = 0; k <= M; k++)
    {
        u = N_f(x);
        if(fabs(u)<std::numeric_limits<double>::epsilon())
        {
            break;
        }
        x = x - u/N_f.diff(x);
    }
    // std::cout << "k = " << k << std::endl;
    // std::cout << "x = " << x << std::endl; 
    return x;
};

class Secant : EquationSolver
{
private:
    double x0;
    double x1;
    double delta;
    int M;
    Function S_f;
public:
    Secant(double _x0, double _x1, double _delta, int _M, Function _S_f);
    double solve();
};

Secant :: Secant(double _x0, double _x1, double _delta, int _M, Function _S_f)
{
    x0 = _x0; x1 = _x1; delta = _delta; M = _M; S_f = _S_f; 
}

double Secant :: solve()
{
    double xn = x1;
    double xn_1 = x0;
    double u = S_f(xn);
    double v = S_f(xn_1);
    int k = 2;
    for( k = 2; k <= M; k++)
    {
        if(fabs(u)>fabs(v))
        {
            double m = xn;
            xn = xn_1;
            xn_1 = m;
            m = u;
            u = v;
            v = m;
        }
        double s = (xn - xn_1)/(u - v);
        xn_1 = xn;
        v = u;
        xn = xn - u*s;
        u = S_f(xn);
        if(fabs(xn - xn_1) < delta || fabs(u) < std::numeric_limits<double>::epsilon())
        {
            break;
        }
    }
    // std::cout << "k = " << k << std::endl;
    // std::cout << "xn_1 = " << xn_1 << std::endl;
    // std::cout << "xn = " << xn << std::endl;
    return xn;
};
